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Abstract 

Timely and accurate information on crop yield is critical to many applications within agriculture 
monitoring. Thanks to its coverage and temporal resolution, coarse spatial resolution satellite 
imagery has always been a source of valuable information for yield forecasting and assessment at 
national and regional scales. With availability of free images acquired by Landsat-8 and Sentinel-2 
remote sensing satellites, it becomes possible to enable temporal resolution of an image every 3—5 
days, and therefore, to develop next generation agriculture products at higher spatial resolution (30 
m). This paper explores the combined use of Landsat-8 and Sentinel-2A for winter crop mapping 
and winter wheat assessment at regional scale. For the former, we adapt a previously developed 
approach for Moderate Resolution Imaging Spectroradiometer (MODIS) at 250 m resolution that 
allows automatic mapping of winter crops taking into account knowledge on crop calendar and 
without ground truth data. For the latter, we use a generalized winter wheat yield model that is 
based on NDVI-peak estimation and MODIS data, and further downscaled to be applicable at 30 m 
resolution. We show that integration of Landsat-8 and Sentinel-2A has a positive impact both for 
winter crop mapping and winter wheat yield assessment. In particular, the error of winter wheat 


yield estimates can be reduced up to 1.8 times comparing to the single satellite usage. 


Keywords: Landsat-8, Sentinel-2, yield, area, mapping, wheat, MODIS, agriculture, Ukraine 


37 


38 


39 


40 


41 


42 


43 


44 


45 


46 


47 


48 


49 


50 


51 


52 


53 


54 


oh) 


56 


57 


58 


59 


60 


61 


62 


1. Introduction 

Timely and accurate information on crop yields at global, national, and regional scales is 
extremely important for many applications [1]. At national/regional scale, it can be an input to local 
authorities to make decisions on food security issues or deciding on subsidies in case of extreme 
weather conditions such as droughts. At field scale, spatial variability of yields can help to obtain 
objective information, for example, for farmers to improve management practices and identify yield 
gaps [2], or for insurance companies to feed this information into insurance models [3, 4]. 

Owing to its coverage, temporal and spatial resolution, remote sensing images from space 
has always been a powerful tool to develop empirical models for predicting and assessing yields at 
regional and national scales [5, 6, 7, 8, 9, 10, 11], or assimilating biophysical parameters into crop 
growth models [12, 13, 14]. In particular, coarse resolution sensors, e.g. Moderate Resolution 
Imaging Spectroradiometer (MODIS), Advanced Very High Resolution Radiometer (AVHRR), 
SPOT-VEGETATION, thanks to its daily coverage and availability of historical datasets, have 
extensively been used for building empirical models for crop yield forecasting and assessment. 
These models connect satellite-derived features, for example vegetation indices (VIs) such as 
Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), Vegetation 
Health Index (VHI) and/or biophysical parameters such as Leaf Area Index (LAI), Fraction of 
Photosynthetically Active Radiation (FPAR), with reference yield data. For example, Johnson 
(2016) [5] analyzed efficiency of multiple MODIS land products including NDVI, EVI, LAI, 
FPAR, and Gross Primary Production (GPP) to assess crop yield at county level in US for ten major 
agriculture commodities. He found positive correlations of vegetation products against yield for all 
crops, except rice, and that finer spatial resolution improved the correlations. Lopez-Lozano et al. 
(2015) [6] investigated the use of the Fraction of Absorbed Photosynthetically Active Radiation 
(fAPAR) derived from SPOT-VEGETATION to assess crop yields (wheat, barley and maize) at 
province level in Europe. They found high correlations (R*>0.6) in water-stressed regions; however, 


lower correlations (R’<0.5) were observed for regions with high yields where water constraints are 


63 


64 


65 


66 


67 


68 


69 


70 


71 


72 


73 


74 


75 


76 


77 


78 


79 


80 


81 


82 


83 


84 


85 


86 


87 


88 


less frequent. Salazar et al. (2007) applied AVHRR-derived VHI to estimate winter wheat yield in 
Kansas, US, and found high correlations with official statistics for 1982—2004 obtaining an error 
around 8%. In order to overcome some limitations of empirical models in terms of robustness, 
Becker-Reshef et al. (2010) [10] developed a generalized winter wheat yield forecasting model that 
was calibrated for one region (US) and successfully applied for another (Ukraine) to provide 
accuracy of less than 10% that is suitable for operational context. Adding meteorological data, in 
particular temperature, has usually had a positive effect on crop yield models reducing the error and 
improving timeliness [5, 6, 7]. Though these models are empirical and based on relative simple 
equations, they perform at the level or even better than more comprehensive crop growth models 
that are based on crop growth simulations [8, 15]. The reasons for that are: complexity of 
accounting multiple factors influencing the yield, lack of high-quality data required to calibrate and 
run such models, and difficulties of upscaling ‘point’ estimates to higher spatial scale [16]. 

Comparing to coarse resolution satellite imagery, the use of Landsat-like (30 m) data to crop 
yield forecasting and assessment has been limited mainly due to lower temporal resolution. 
Nevertheless, there were studies aiming at fusing Landsat with MODIS data [17, 18], and 
combining Landsat with biophysical models [19, 20]. However, these approaches showed mixed 
results in terms of errors and still had limitations constrained by lower frequency of moderate 
resolution images. With the combined use of Landsat-8 and Sentinel-2 remote sensing satellites that 
would enable acquisition of an image every 3—5 days globally, as well as development of cloud 
platforms such as Google Earth Engine (GEE) [20, 21, 22], it becomes possible to implement 
approaches similar to those used for MODIS/AVHRR to develop next generation agriculture 
products at higher spatial resolution (30 m). 

This paper presents one of the first studies to combine Landsat-8 and Sentinel-2A imagery 
for crop yield mapping by downscaling a generalized empirical model developed for MODIS data 
[7, 10]. The model is based on capturing the peak NDVI to correlate with the yield, and growing 


degree days (GDD) to improve the timeliness of the model. Therefore, the main objectives of the 


89 


90 


91 


92 


93 


94 


95 


96 


97 


98 


99 


100 


101 


102 


103 


104 


105 


106 


107 


108 


109 


110 


111 


112 


113 


study are: (i) to assess performance of downscaling a generalized NDVI-based empirical model for 
winter wheat yield forecasting from coarse spatial resolution to moderate one at 30 m; (ii) to 
explore the combined use of images acquired by Landsat-8 and Sentinel-2A remote sensing 


satellites for winter crop mapping and winter wheat yield assessment at regional level. 


2. Study area & materials 
2.1. Study area and reference data 

The study is performed for Kirohohradska oblast in Ukraine for 2016 (Fig. 1). Oblast is a 
high-level administrative division of the country (there are 24 oblasts in Ukraine and Autonomous 
Republic of Crimea), and each oblast is further divided into districts. Kirovhradska oblast is located 
in the central part of Ukraine and composed of 21 districts with geographical area ranging from 65 
to 165 thousand ha and cropland area ranging from 27 to 112 thousand ha. The reasons for selecting 
this region is that it is a top 10 wheat producer in Ukraine and because of availability of reference 
crop yield and harvested area data at district scale. Winter wheat is one of the major crops in 
Kirovhradska oblast accounting for 20% of production of all crops in the region. Winter wheat is 
mainly rain-fed in the region and usually planted in September-October. After dormancy during the 
winter, it emerges early spring reaching maturity by the end of June. Harvest of winter crops is 
typically undertaken in July. 

Reference data on crop yield and harvested area at district level were collected from the 
Department of Agro-Industry Development of Kirovohrad State Administration (http://apk.kr- 
admin.gov.ua). The data were made available online as the harvest progressed and were based on 
farm surveys of all large agricultural enterprises (that account of more than 90% of all winter crops 
production in the region) and samples of household farms the same way as official statistics is 
collected [23]. The final estimates for winter crop yields and areas were available at the end of 


November and were used as reference in this study. 
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Fig. 1. A map of Ukraine with division into administrative regions (oblasts). The study area 
(Kirovohradska oblast) is highlighted on /eft figure and shown with division into district on right 


figure. 


2.2. Landsat-8/OLI and Sentinel-2A /MSI datasets 

Remote sensing images acquired by the Operational Land Imager (OLI) instrument aboard 
of Landsat-8 satellite and by the Multi-Spectral Instrument (MSI) aboard Sentinel-2A satellite were 
used in the study. Landsat-8/OLI captures images of the Earth’s surface in 9 spectral bands at 30 m 
spatial resolution (15 m for panchromatic band) [24] while Sentinel-2A/MSI captures images of the 
Earth’s surface in 13 spectral bands at 10 m, 20 m and 60 m spatial resolution [25]. Overall, 51 
Landsat-8 and 87 Sentinel-2A scenes were acquired over the study area from March 1, 2016 to July, 
31, 2016. Landsat-8 scenes covered the following coordinates (path/row) of the World-wide 
Reference System (WRS-2): 178/026, 179/026, 179/027, 180/026, 180/027, and 181/026. The 
swath of the Landsat-8 scene is approximately 185 km x 180 km. Sentinel-2A scenes covered the 
following tiles: 35UQQ, 35UQP, 36UUV, 36UUU, 36UVV, 36UVU, 36UWV, and 36UWU. The 
size of the Sentinle-2A tile is approximately 110 km x 110 km. 

The Landsat-8/OLI and Sentinel-2A/MSI scenes were atmospherically corrected for surface 
reflectance using the LaSRC algorithm [26] (Fig. 2) ensuring consistency between these datasets as 
well as with MODIS data used for building a generalized crop yield model [10, 28]. Cloud and 
shadow screening for Landsat-8 and Sentinel-2A scenes was performed using the Fmask algorithm 


[27] and inversion residuals from aerosol optical thickness (AOT) estimation [26] (Fig. 3). The 
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pixels identified as those with high aerosol content were also masked out. Images from Sentinel- 
2A/MSI were further converted to 30 m to match spatial resolution of Landsat-8/OLI. Since 
atmospheric correction for Sentinel-2A was performed at 10 m spatial resolution for all spectral 
bands, conversion to 30 m was carried out by aggregation (averaging). 

It was found that Landsat-8/OLI and Sentinel-2A/MSI exhibit misregistration issues [29]; 
therefore additional co-registration was performed to ensure spatial consistency between the 
datasets [30]. Finally, NDVI was calculated for Landsat-8 scenes using band 5 (near-infra red — 


NIR) and band 4 (red), and for Sentinel-2A scenes using band 8A (NIR) and band 4 (red). 


Landsat-8, 30 m, true color (4-3-2), 17-Apr-2016 Sentinel-2A, 30 m, true color (4-3-2), 18-Apr-2016 


Landsat-8, 30 m, false color (5-4-3), 17-Apr-2016 Sentinel-2A, 30 m, false color (8a-4-3), 18-Apr-2016 
" 7 ‘ 


Fig. 2. Examples of images acquired by Landsat-8 and Sentinel-2A satellites 1 day apart and 
atmospherically corrected using the LaSRC algorithm [26]. True colour images were composed of 
bands 4-3-2 for Landsat-8 and Sentinel-2A, and scaled from 0 to 0.15. False colour images were 
composed of bands 5-4-3 for Landsat-8 and 8A-4-3 for Sentinel-2A, and scaled from 0 to 0.3 for 


NIR, and 0 to 0.1 for red and green bands. 


Sentinel-2A, 30 m, true color Sentinel-2A, high residual 
(4-3-2), 17-Jun-2016 (shown in black) 
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151 
152. Fig. 3. Example of cloud and shadow detection for Sentinel-2A images 

153 

154 3. Methodology 

155 Winter wheat yield mapping and assessment at regional scale consists of the two major 
156 steps: (1) winter crop mapping; (ii) yield assessment at 30 m spatial resolution. Fig. 4 illustrates all 
157. processing steps along with input datasets. These steps are described in detail in the following sub- 


158 sections. 


Landsat-8 
Sentinel-2A 


Calculate maximum Calculate maximum 
NDVI from March 1 NDVI from March 1 
to April 6 to June 30 


Cropland pana Run yield model 
anask model over cropland with slope 4.61 
mask (Eq. 1) . 
Winter crop 
map 


160 Fig. 4. Algorithm flowchart. 
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3.1. Winter crop mapping 

For winter crop mapping, we adopted a previously developed approach for MODIS [31] that 
allows automatic mapping of winter crops using a priori knowledge on crop calendar and without 
using reference (ground truth) data. The method is based on per-pixel estimation of the peak NDVI 
(metric) during early spring (or early fall depending on hemisphere), when winter crops have 
developed biomass, while other crops (spring and summer) have no biomass in that time period. 
The calculated metric will have high NDVI values for winter crops and low NDVI values for other 
crops (Fig. 5). Then, the metric is modelled using a Gaussian mixture model (GMM) [32] to 
automatically discriminate different crop types (winter versus others). The GMM is a linear 
combination of Gaussian distributions that can model any continuous distribution: 
p(X) = Ditar TN loa, Ded (1) 
where each Gaussian density N(x|uU,, ¥,) is called a component of the mixture and has its own 


mean pl, and covariance Y},; parameters 7, are weight (mixing) coefficients with )4_,2, =1 
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Fig. 5. Empirical distribution for the peak NDVI and fitted GMM model. The solid green line 


shows the fitted GMM distribution, while the dashed lines show the mixture model components. 
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Parameters of the GMM model are estimated using an expectation-maximization (EM) 
algorithm that is run for all pixels identified as cropland. In our study, we used a cropland layer 
acquired from the land cover map generated for Ukraine at 30 m spatial resolution [33]. The 
constraint to utilize cropland pixels only comes from potential confusion with grassland, hay, 
bulrush that might also have already developed biomass within indicated time period. The 
component with the largest mean in the obtained GMM model is considered to belong to the winter 
crop class (Fig. 5). Finally, the derived GMM model is applied to all cropland pixels, and a 
posteriori probability (Eq. 1) of the pixel belonging to the winter crop class is estimated in the final 


resulting map. Pixels with probability larger than 0.5 are considered as winter crops. 


3.2. Winter wheat yield mapping and assessment 

Peak NDVI estimated on a per-pixel basis from Landsat-8/OLI and Sentinel-2A/MSI images 
from March to June was selected as a primary parameter for assessing winter wheat yield. In 
multiple studies NDVI has been shown to be strongly correlated with yields for a variety of crop 
types [5, 8, 9, 10]. Since there are no available historical data for a combination of Landsat-8 and 
Sentinel-2A images to correlate with yield measurements and build a crop yield model at district 
scale, we used a MODIS-derived winter wheat yield model that was calibrated for US and directly 
applied for Ukraine [7, 10]. More specifically, the model takes advantage of daily MODIS data at 
Climate Modeling Grid (CMG) scale at 0.05° resolution to capture an NDVI peak and correlate 
with the yield. However, since proportion of winter wheat is variable within the CMG pixel, the 
model establishes a generalized relationship between the slope of NDVI against yield and pixel 
purity [10]: s=9.61—0.05*m, where m is the winter wheat proportion at CMG scale (from 0 to 
100%), and s is the slope such as yield=s*NDVI. 

In case of Landsat-8—Sentinel-2A images, we can assume that purity at 30 m level is 100%, 
i.e. m=100. Therefore, we obtain the slope of 4.61 to be applied to an NDVI peak calculated from 


the combination of Landsat-8 and Sentinel-2A data to map winter wheat yield at 30 m resolution. 
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Therefore, the MODIS-derived coarse resolution (0.05°) winter wheat yield model, that was 
calibrated for US [10], is downscaled using winter wheat purity as a proxy to derive the slope 
between the peak NDVI and yield. This slope (4.61) is directly applied to the peak NDVI calculated 
from Landsat-8—Sentinel-2A images to derive a winter wheat yield map at 30 m resolution. These 
are used to estimate district-level yields by averaging yields at 30 m over winter crop masks 
(section 3.1) for each district. In addition to the average, a standard deviation and coefficient of 
variation (CV), defined as a ratio between the standard deviation and the mean, is estimated as well. 
The estimated district-level yields are validated using independent reference data (section 2.1) 


collected at district level in Kirvohradska oblast in Ukraine. 


3.3. Validation metrics 
For comparison of satellite-derived winter crop areas and winter wheat yield with reference 
datasets at district level, we used the APU analysis metrics [28]: 


* accuracy (A) that shows the average bias of the estimates 
1 
A= ei = 0;), (2) 


* precision (P) that shows repeatability of the estimates 


(POH ADs, (3) 


* uncertainty (U) that is the root mean squared error 


U= [EON - 0%, (4) 


¢ relative uncertainty (rU) normalized by an average of reference values: 


rU(%) =7~— x100%, (5) 
yLi=1 0; 


where P; and O; are computed (from satellites) and observed (from reference) values, respectively. 
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231 4. Results & discussion 

232. = 4.1. Winter crop mapping 

233 The GMM approach to winter crop mapping was applied to the peak NDVI calculated for 
234 the time period from March | to April 6 using a combination of Landsat-8 and Sentinel-2A, as well 
235 as using each of them separately. This was done in order to assess an added value of the combined 
236 use of these datasets. The indicated period (March | to April 6) was selected in such a way to 
237. ~capture NDVI development of winter crops and avoid confusion with early spring cereals that were 
238 planted beginning of March. The derived maps were used to calculate the area of winter crops at 
239 districts level by pixel-counting. These estimates were compared to reference values and are 
240 presented in Table 1 and Fig. 6. The derived winter crop using Landsat-8 and Sentinel-2A is 
241 ~— illustrated in Fig. 7. 

242 

243 Table 1. Comparison of satellite-derived winter crop areas with official statistics on harvested areas 


244 at district level. Estimates of the APU metrics are given in ha. 


Metric LC8-S2A LC8 S2A 
A 612 1081 839 
P 1719 5061 1962 
U 1785 5056 2090 
rU, % 11.6 32.7 13.5 
R 0.90 0.64 0.88 
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Fig. 6. Plots of official statistics on harvested winter crop areas against satellite-derived ones using 


a combination of Landsat-8 and Sentinel-2A (/eft), Landsat-8 only (centre), and Sentinel-2A only 


(right). 
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Fig. 7. Final map of winter crops derived from Landsat-8 and Sentinel-2A images using the GMM 


approach for Kirohradska oblast in 2016. 


Combination of Landsat-8 and Sentinel-2A allowed us to achieve R’=0.9 and relative 
uncertainty of 11.6% when estimating winter crop areas at district level. It should be noted that 
these results were achieved in an automatic way utilizing knowledge on crop calendar and without 
utilizing any ground truth data. The use of Landsat-8 images only did not produce satisfactory 
results (R’ = 0.64 and relative uncertainty of 32.7%) because of unavailability of cloud-free images 


early spring especially in the eastern districts of the oblast whereas the use of Sentinel-2A yielded 
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263 R*=0.88 and relative uncertainty of 13.5%. Overall, these results demonstrate the benefits, in a 
264 quantitative way, of the combined use of Landsat-8 and Sentinel-2A satellites comparing to the 
265 single-satellite usage. 

266 

267 4.2. Winter wheat yield mapping 

268 Results of comparison of the estimated winter wheat yields at district level are presented in 
269 Table 2 and Fig. 8. 

270 

271. Table 2. Comparison of satellite-dertved winter wheat yields with official statistics at district level. 


272 Estimates of the APU metrics are given in t/ha. 


Metric LC8-S2A LC8 S2A 

A -0.17 -0.48 -0.34 

Pr 0.26 0.31 0.32 

U 0.31 0.57 0.46 

rU, % TA 14.3 11.5 

R 0.45 0.29 0.28 
273 
274 as tad 5 ce mha ” si re . >> sat, rina“? = . sat, na *° 
275 ~=~*Fig. 8. Plots of official statistics on winter wheat yield against satellite-derived ones using a 
276 combination of Landsat-8 and Sentinel-2A (left), Landsat-8 only (centre), and Sentinel-2A only 
277 ~—s (right). 
278 
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As with winter crop areas, the combination of Landsat-8 and Sentinel-2A outperformed the 
single satellite usage. When using either Landsat-8 or Sentinel-2A, the peak NDVI approach 
underestimated official statistics by -0.48 t/ha and -0.34 t/ha, respectively, while their combination 
improved to -0.17 t/ha. In terms of uncertainty, the peak NDVI approach for the Landsat-8— 
Sentinel-2A combination provided 0.31 t/ha (7.7%) whereas those values were 1.8 times higher for 
the Landsat-8 usage only (0.57 t/ha, 14.3%) and 1.5 times higher for the Sentinel-2A usage only 
(0.46 t/ha, 11.5%). These results clearly demonstrate the importance of higher observation 
frequency achieved with combination of Landsat-8 and Sentinel-2A satellites comparing to the 
single use. 

The results presented in Fig. 8 (left) were further analyzed for errors. Overall, the points can 
be divided into 3 groups. The first group is composed of 3 points representing districts with official 
statistics yields values close to 4 t/ha and underestimated by the peak NDVI approach. These 
districts feature relatively large values of CV of 21% whereas the average CV for all other districts 
is approximately 13%. The reason for that is smaller number of images available for these districts 
(mainly in the eastern part) which reduces ability to capture the peak NDVI. The second group is 
composed of districts with official statistics yields larger than 4 t/ha. The reason for that is 
saturation of NDVI occurs and the proposed approach fails to discriminate yield values at this level. 
Fig. 9 shows an example of NDVI time-series for the district with reference yield of 4.3 t/ha and 
estimated yield of 4.04+0.40 t/ha with NDVI quickly achieving the value of 0.8 on April 29 (day of 
the year (DOY) 120) and not changing considerably (within 0.8—0.9) during the following days 50 
days (until June 18 or DOY=170). The NDVI values start to decrease when the senescence phase 


occurs and the crop is eventually harvested. 
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Fig. 9. A combined Landsat-8—Sentinel-2A derived NDVI time-series of winter wheat for the 


district with reference yield at 4.3 t/ha. 


The third group involves 8 districts with moderate yield values of up to 4 t/ha. The proposed 
approach is able to explain variations in the winter wheat yield (R’=0.8) giving a bias of 0.1 t/ha 


and uncertainty of U=0.13 t/ha (3.5%). 


6. Conclusion 

This study attempted to explore the combined use of Landsat-8 and Sentinel-2A satellites to 
winter crop mapping and winter wheat yield assessment at regional level. For both tasks, the 
increased frequency of observations from the satellites was critical as it allowed us to achieve better 
performance comparing to the single satellite usage. For winter crop mapping, we adopted a 
previously developed approach for MODIS that allowed automatic winter crop mapping taking into 
account a priori knowledge on crop calendar without utilizing ground reference data. When 
comparing to official statistics on winter crop harvested areas, this approach gave R°=0.9 and 
relative error of 11.6%. These results are encouraging as with little data inputs (crop calendar and 
cropland mask) and high temporal resolution of Landsat-8—Sentinel-2A satellites, it would allow the 


creation of winter crop maps at global scale at 30 m resolution. 
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For winter wheat yield mapping, we downscaled the generalized empirical model that is 
based on peak NDVI and was developed using MODIS data, and directly applied it to the Landsat- 
8—Sentinel-2A images. The model was efficient in explaining moderate yield values (<4 t/ha) with 
R?=0.8; however, it failed to capture the variance of high yield values (>4 t/ha) due to NDVI 
saturation. Overall, the downscaled peak NDVI approach with combined use of Landsat-8 and 
Sentinel-2A gave uncertainty of 0.31 t/ha (7.7%) and R°=0.45 substantially outperforming Landsat- 


8 only (1.8 times) and Sentinel-2A only (1.5 times). 
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